Influence of Composite Structure on Temperature Distribution—An Analysis Using the Finite Difference Method

Among composites, we can distinguish periodic structures, biperiodic structures, and structures with a functional gradation of material properties made of two or more materials. The selection of the composite’s constituent materials and the way they are distributed affects the weight of the composite, its strength, and other properties, as well as the way it conducts heat. This work is about studying the temperature distribution in composites, depending on the type of component material and its location. For this purpose, the Tolerance Averaging Technique and the Finite Difference Method were used. Differential equations describing heat conduction phenomena were obtained using the Tolerance Averaging Technique, while the Finite Difference Method was used to solve them. In terms of results, temperature distribution plots were produced showing the effect of the structure of the composite on the heat transfer properties.


Introduction
Composite materials are materials with specific designed properties. These properties can be obtained by combining certain materials in the relevant way. Composite structures can be characterized by better strength, resistance, or thermal conductivity than their individual component materials. In periodic and biperiodic structures or in those with a functional gradation of properties, it is possible to separate mentally repeatable cells with the same structure. The dimension of such a cell is called the microstructure parameter. Methods that are commonly used in such considerations are certainly the Finite Element Method [1,2]; asymptotic homogenization [3]; the variation of homogenization, introducing additional unknowns called microlocal parameters [4,5]; the higher order theory [6,7]; and the boundary element method [8]. Not all of them take into account the microstructure parameter, in contrast to the Tolerance Averaging Technique, which is used in this work [9,10]. This method is used in the analysis of a variety of structures, such as beams [11,12], plates [13] and shells [14,15]. The main reason for utilizing the Tolerance Averaging Technique in this study was the great alignment of its results to the exact solution when it comes to the stationary problem of heat conduction [16].
The Tolerance Averaging Technique makes it possible to average out the discontinuous coefficients in the heat conduction equation, which result from the different properties of the composite's constituent materials. The resulting averaged equations, describing the heat flow in the composite structures, were solved numerically using the Finite Difference Method, which was implemented by the authors into MAPLE 2019 software. The proposed FDM algorithm is universal in terms of the dimensions of the composite, its structure, and the type of constituent materials, but not in terms of the assumed boundary conditions of thermal loads, because the number of equations of the Finite Difference Method depends on the number of nodes wherein the sought unknown function is not predetermined. predetermined.

Various Composite Structures of Interest
Composite materials are often used in the design of plate and shell structures. different structures of this type were analyzed. The first structure is a periodic com as shown in Figure 1a. The dimensions of the composite are denoted by L1 and L2 a material properties change recurrently along the x2-direction.
The dimension of the mentally separated cell, made of two materials, is deno l2. The volume proportion of individual materials in the cell is constant. The volum portion of the first material is denoted by v1, and the volume share of the second m is denoted by v2, while v1 + v2 = 1.
The second structure is a structure with a functional gradation of properties shown in Figure 1b. The volume proportion of individual materials in the cell is not constant and termined by a function that depends on the x2-coordinate. Therefore, the volume p tion of the first material is denoted by v1(x2), and the volume share of the second m is denoted by v2(x2).
The last one-the third composite-is a biperiodic structure (refer to Figure 2). case of this structure, the properties change recurrently along both directions-x1 a The dimensions of the mentally separated cell are denoted by l1 and l2.  The dimension of the mentally separated cell, made of two materials, is denoted by l 2 . The volume proportion of individual materials in the cell is constant. The volume proportion of the first material is denoted by v 1, and the volume share of the second material is denoted by v 2 , while v 1 + v 2 = 1.
The second structure is a structure with a functional gradation of properties and is shown in Figure 1b.
The volume proportion of individual materials in the cell is not constant and is determined by a function that depends on the x 2 -coordinate. Therefore, the volume proportion of the first material is denoted by v 1 (x 2 ), and the volume share of the second material is denoted by v 2 (x 2 ).
The last one-the third composite-is a biperiodic structure (refer to Figure 2). In the case of this structure, the properties change recurrently along both directions-x 1 and x 2 . The dimensions of the mentally separated cell are denoted by l 1 and l 2 . Method depends on the number of nodes wherein the sought unknown function is not predetermined.

Various Composite Structures of Interest
Composite materials are often used in the design of plate and shell structures. Three different structures of this type were analyzed. The first structure is a periodic composite, as shown in Figure 1a. The dimensions of the composite are denoted by L1 and L2 and the material properties change recurrently along the x2-direction.
The dimension of the mentally separated cell, made of two materials, is denoted by l2. The volume proportion of individual materials in the cell is constant. The volume proportion of the first material is denoted by v1, and the volume share of the second material is denoted by v2, while v1 + v2 = 1.
The second structure is a structure with a functional gradation of properties and is shown in Figure 1b. The volume proportion of individual materials in the cell is not constant and is determined by a function that depends on the x2-coordinate. Therefore, the volume proportion of the first material is denoted by v1(x2), and the volume share of the second material is denoted by v2(x2).
The last one-the third composite-is a biperiodic structure (refer to Figure 2). In the case of this structure, the properties change recurrently along both directions-x1 and x2. The dimensions of the mentally separated cell are denoted by l1 and l2.  The dimension of the cell, called the microstructure parameter, is dependent on the number of composite cells N.

Tolerance Averaging Technique
The Tolerance Averaging Technique, mentioned previously, was used to average Equation (1) of the heat transfer issue [17] to replace the terms with discontinuous coefficients: where θ symbolizes the total temperature field; the components of the tensor of conductivity are denoted by k ij (i,j = 1, 2, 3), a specific heat by c and a mass density by ρ. Each cell is composed of two (periodicity) or four (biperiodicity) sub-cells, and the above-mentioned material properties within a given sub-cell are constant. The technique, used to average the discontinuous coefficients in Equation (1), is constantly developed and used in the analysis of periodic [18][19][20][21][22][23], biperiodic [24][25][26] and functionally graded structures [27,28]. The most relevant assumption of the Tolerance Averaging Technique, from the point of view of the carried-out considerations, is the micro-macro decomposition assumption, according to Equation (2): where the total temperature field θ is divided into an averaged part ϑ and an oscillating part in the form of product g a ·ψ a . The average temperature (the so-called macrotemperature) is denoted by ϑ, the fluctuation amplitude of the temperature (the new additional unknowns) is denoted by ψ a , and the fluctuation shape functions are denoted by g a . With reference to periodic structures and structures with a functional gradation of properties, index a in (2) is equal to one, and the saw-type function is assumed as the g 1 -function. The g 1 -function was chosen based on the available literature and taking into account the specifics of the issue under consideration (the heat conduction issue). This function, shown in Figure 3, provides the continuity of the temperature field and the heat flux between the cells and on interfaces between the first and the second material. The dimension of the cell, called the microstructure parameter, is dependent on the number of composite cells N.

Tolerance Averaging Technique
The Tolerance Averaging Technique, mentioned previously, was used to average Equation (1) of the heat transfer issue [17] to replace the terms with discontinuous coefficients: where θ symbolizes the total temperature field; the components of the tensor of conductivity are denoted by kij (i,j = 1, 2, 3), a specific heat by c and a mass density by ρ. Each cell is composed of two (periodicity) or four (biperiodicity) sub-cells, and the above-mentioned material properties within a given sub-cell are constant. The technique, used to average the discontinuous coefficients in Equation (1), is constantly developed and used in the analysis of periodic [18][19][20][21][22][23], biperiodic [24][25][26] and functionally graded structures [27,28]. The most relevant assumption of the Tolerance Averaging Technique, from the point of view of the carried-out considerations, is the micromacro decomposition assumption, according to Equation (2): where the total temperature field θ is divided into an averaged part ϑ and an oscillating part in the form of product ga·ψa. The average temperature (the so-called macrotemperature) is denoted by ϑ, the fluctuation amplitude of the temperature (the new additional unknowns) is denoted by ψa, and the fluctuation shape functions are denoted by ga. With reference to periodic structures and structures with a functional gradation of properties, index a in (2) is equal to one, and the saw-type function is assumed as the g1-function. The g1-function was chosen based on the available literature and taking into account the specifics of the issue under consideration (the heat conduction issue). This function, shown in Figure 3, provides the continuity of the temperature field and the heat flux between the cells and on interfaces between the first and the second material. With reference to biperiodic structure index, a in (2) runs over 1 and 2. For this structure, which has the fluctuation shape functions g1 and g2, some combination of the sawtype fb and piecewise parabolic hb function was chosen, and it is shown in Figure 4. The fluctuation shape function g1(x1,x2) can be expressed as a product of functions f1(x1) and With reference to biperiodic structure index, a in (2) runs over 1 and 2. For this structure, which has the fluctuation shape functions g 1 and g 2 , some combination of the saw-type f b and piecewise parabolic h b function was chosen, and it is shown in Figure 4. The fluctuation shape function g 1 (x 1 ,x 2 ) can be expressed as a product of functions f 1 (x 1 ) and f 2 (x 2 ), and the fluctuation shape function g 2 (x 1 ,x 2 ) can be expressed as a product of functions h 1 (x 1 ) and h 2 (x 2 ). f2(x2), and the fluctuation shape function g2(x1,x2) can be expressed as a product of functions h1(x1) and h2(x2). In addition to the micro-macro decomposition assumption, the Tolerance Averaging Technique also introduces the other assumptions and definitions such as the averaging operation, tolerance-periodic function, slowly varying function (examples include the macro-temperature function ϑ and the fluctuation amplitude ψi) and the highly oscillating function (an example is the fluctuation shape function ga).
By using the above definitions and assumptions of the Tolerance Averaging Technique discussed in [29], the equations with the averaged coefficients describing the heat transfer phenomenon were obtained. The equations in reference to the periodic structure and the structure with a functional gradation of properties (the so-called tolerance model equations) are shown below: where tensor, whose components are kij, is denoted by K, and is the tensor of conductivity; gradient operator ∇ is expressed as follows (∂1,∂2,∂3); operator ∂ = (∂1,0,0) is a gradient in the x1-direction; and the overlined nabla operator is a gradient in the x2-and x3-directions. The equations in reference to the biperiodic structure are shown below:  In addition to the micro-macro decomposition assumption, the Tolerance Averaging Technique also introduces the other assumptions and definitions such as the averaging operation, tolerance-periodic function, slowly varying function (examples include the macro-temperature function ϑ and the fluctuation amplitude ψ i ) and the highly oscillating function (an example is the fluctuation shape function g a ).
By using the above definitions and assumptions of the Tolerance Averaging Technique discussed in [29], the equations with the averaged coefficients describing the heat transfer phenomenon were obtained. The equations in reference to the periodic structure and the structure with a functional gradation of properties (the so-called tolerance model equations) are shown below: where tensor, whose components are k ij , is denoted by K, and is the tensor of conductivity; gradient operator ∇ is expressed as follows (∂ 1 ,∂ 2 ,∂ 3 ); operator ∂ = (∂ 1 ,0,0) is a gradient in the x 1 -direction; and the overlined nabla operator is a gradient in the x 2 -and x 3 -directions. The equations in reference to the biperiodic structure are shown below: where operator ∂ = (∂ 1 ,∂ 2 ,0) is a gradient in the x 1 -and x 2 -directions and the overlined ∇ operator is a gradient in the x 3 -direction.

Finite Difference Method
It is well-known that the Finite Difference Method, when applied to the differential equation (DE), ordinary equation (ODE) or partial equation (PDE), leads to the system of algebraic linear equations, which is its biggest advantage. The domain of unknown function is narrowed down to a grid of nodes, while the codomain is replaced by set of function values in these nodes. Moreover, every derivative appearing in the DE is approximated by the finite difference quotient, depending very often on function values at adjacent nodes, and that, in particular, may cause a need to temporarily create and use virtual nodes (outside the domain). Function values at virtual nodes are acquired from boundary conditions.
What distinguishes the considered problem from the others is that we actually deal here with the system of PDEs, thus we have more than one unknown function to find. Moreover, each one may have different type of condition in common for all boundaries. Hence, most of the work that has been conducted here was grouping and ordering unknowns and constructing matrix of coefficients for the system. When everything was settled, the MAPLE environment were used in order to implement the Finite Difference Method, and that implementation is further called an algorithm.
To create an algorithm of the Finite Difference Method, it is necessary to divide an analyzed structure into ranges, redefine the tolerance model equations to formulate them in index notation, replace the derivatives appearing in these equations by appropriate differential quotients, and define the boundary conditions, because this affects the number of equations to solve [26]. All of the structures, whether periodic, biperiodic, or with a functional gradation of properties, were discretized along the x 1 -direction into m ranges, so the number of the nodes along the x 1 -direction equals m + 1, and analogously along the x 2 -direction into n ranges, so the number of the nodes along the x 2 -direction equals n + 1; refer to Figure 5.  ρ ϑ + ρ ψ + ρ ψ + ⋅∂ ⋅∇ϑ+ ⋅∂ ∂ ψ + ⋅∂ ⋅∇ψ c g g c g g g g g g g g gg gg g g g g g g g g where operator ∂ = (∂1,∂2,0) is a gradient in the x1-and x2-directions and the overlined ∇ operator is a gradient in the x3-direction.

Finite Difference Method
It is well-known that the Finite Difference Method, when applied to the differential equation (DE), ordinary equation (ODE) or partial equation (PDE), leads to the system of algebraic linear equations, which is its biggest advantage. The domain of unknown function is narrowed down to a grid of nodes, while the codomain is replaced by set of function values in these nodes. Moreover, every derivative appearing in the DE is approximated by the finite difference quotient, depending very often on function values at adjacent nodes, and that, in particular, may cause a need to temporarily create and use virtual nodes (outside the domain). Function values at virtual nodes are acquired from boundary conditions.
What distinguishes the considered problem from the others is that we actually deal here with the system of PDEs, thus we have more than one unknown function to find. Moreover, each one may have different type of condition in common for all boundaries. Hence, most of the work that has been conducted here was grouping and ordering unknowns and constructing matrix of coefficients for the system. When everything was settled, the MAPLE environment were used in order to implement the Finite Difference Method, and that implementation is further called an algorithm.
To create an algorithm of the Finite Difference Method, it is necessary to divide an analyzed structure into ranges, redefine the tolerance model equations to formulate them in index notation, replace the derivatives appearing in these equations by appropriate differential quotients, and define the boundary conditions, because this affects the number of equations to solve [26]. All of the structures, whether periodic, biperiodic, or with a functional gradation of properties, were discretized along the x1-direction into m ranges, so the number of the nodes along the x1-direction equals m + 1, and analogously along the x2-direction into n ranges, so the number of the nodes along the x2-direction equals n + 1; refer to Figure 5. Then, the thermal boundary conditions were formulated, dividing the composite into the areas shown in Figure 6. Thus, ∆x 1 equals L 1 /m and ∆x 2 equals L 2 /n, where L 1 and L 2 are the dimensions of the composite.
Then, the thermal boundary conditions were formulated, dividing the composite into the areas shown in Figure 6.
At the nodes of area 1 (i equals 1 and j is from 1 to n + 1), the total temperature θ is known (defined) and hence the macro-temperature ϑ is also known at these nodes. Analogously, at the nodes of area 2 (i is from 2 to m + 1 and j equals 1), the macrotemperature ϑ is known. The nodes of area 3 (i is from 2 to m and j equals n + 1) and the nodes of area 4 (i equals m + 1 and j is from 2 to n) belong to the right and bottom edges of the composite, which were assumed as thermally insulated, which leads to the following conditions: At the nodes of area 1 (i equals 1 and j is from 1 to n + 1), the total temperature θ is known (defined) and hence the macro-temperature ϑ is also known at these nodes. Analogously, at the nodes of area 2 (i is from 2 to m + 1 and j equals 1), the macro-temperature ϑ is known. The nodes of area 3 (i is from 2 to m and j equals n + 1) and the nodes of area 4 (i equals m + 1 and j is from 2 to n) belong to the right and bottom edges of the composite, which were assumed as thermally insulated, which leads to the following conditions: At the nodes of area 5 (i equals m + 1 and j equals n + 1), both conditions expressed by Equations (8)-(9) are defined. Then, the boundary conditions for fluctuation amplitudes of the temperature ψ1 and ψ2 are defined, following the areas shown in Figure 7. At the nodes of area 1 (i equals 1 and j is from 1 to n + 1), the fluctuation amplitudes of the temperature ψ1 and fluctuations amplitudes of the temperature ψ2 are defined. Likewise, the conditions in other areas are: area 2 (i is from 2 to m and j equals 1), area 3 (i is from 2 to m and j equals n + 1) and area 4 (i equals m + 1 and j is from 1 to n + 1). After assuming the boundary conditions, the matrix of coefficients appearing in the equations of the tolerance model was defined. These coefficients can be grouped in At the nodes of area 5 (i equals m + 1 and j equals n + 1), both conditions expressed by Equations (8)-(9) are defined. Then, the boundary conditions for fluctuation amplitudes of the temperature ψ 1 and ψ 2 are defined, following the areas shown in Figure 7. At the nodes of area 1 (i equals 1 and j is from 1 to n + 1), the fluctuation amplitudes of the temperature ψ 1 and fluctuations amplitudes of the temperature ψ 2 are defined. Likewise, the conditions in other areas are: area 2 (i is from 2 to m and j equals 1), area 3 (i is from 2 to m and j equals n + 1) and area 4 (i equals m + 1 and j is from 1 to n + 1). At the nodes of area 1 (i equals 1 and j is from 1 to n + 1), the total temperature θ is known (defined) and hence the macro-temperature ϑ is also known at these nodes. Analogously, at the nodes of area 2 (i is from 2 to m + 1 and j equals 1), the macro-temperature ϑ is known. The nodes of area 3 (i is from 2 to m and j equals n + 1) and the nodes of area 4 (i equals m + 1 and j is from 2 to n) belong to the right and bottom edges of the composite, which were assumed as thermally insulated, which leads to the following conditions: At the nodes of area 5 (i equals m + 1 and j equals n + 1), both conditions expressed by Equations (8)-(9) are defined. Then, the boundary conditions for fluctuation amplitudes of the temperature ψ1 and ψ2 are defined, following the areas shown in Figure 7. At the nodes of area 1 (i equals 1 and j is from 1 to n + 1), the fluctuation amplitudes of the temperature ψ1 and fluctuations amplitudes of the temperature ψ2 are defined. Likewise, the conditions in other areas are: area 2 (i is from 2 to m and j equals 1), area 3 (i is from 2 to m and j equals n + 1) and area 4 (i equals m + 1 and j is from 1 to n + 1). After assuming the boundary conditions, the matrix of coefficients appearing in the equations of the tolerance model was defined. These coefficients can be grouped in After assuming the boundary conditions, the matrix of coefficients appearing in the equations of the tolerance model was defined. These coefficients can be grouped in reference to the equations and individual node areas [25,26]. In an analogous way, the vector of free terms is defined. These terms result from the defined boundary conditions, which are the known values of the macro-temperature ϑ on the top and the left edges of the composite (area 1 and area 2; refer to Figure 6) and the known values of the fluctuation amplitudes of the temperature ψ 1 and ψ 2 on all of the edges of the composite (areas 1 to 4; refer to Figure 7). Thus, a system of non-uniform equations was formulated. To solve this system, some variation of the Finite Difference Method was used-the so-called Crank-Nicolson method (the mixed method). The characteristic parameter for this method was assumed to be equal to a half. The Crank-Nicolson method guarantees the convergence of the solution regardless of the division of the structure.

Results
Using the created algorithm of the Finite Difference Method [23,25,26,28], calculations were conducted for the three structures described above. The two-dimensional, non-stationary heat-conduction problem of composites with base dimensions equal to L 1 = L 2 = 1 [m] was analyzed. During the analysis, the composite was assumed to consist of ten cells. In the case of the periodic structure, the dimension of the first sub-cell γ 2 is taken initially as equal to 0.25; however, in the case of the structure with a functional gradation of properties, the function of gradation γ 2 (x 2 ) equals x 2 /L 2 . On the other hand, in reference to the biperiodic structure, the volume share of the first sub-cell is originally equal to γ 1 · γ 2 , where γ 1 = γ 2 = 0. 25

Periodic and Biperiodic Structure
Firstly, in this section, the approximated results of the Finite Difference Method are presented in the form of 3D-maps of the averaged temperature ϑ and the fluctuation amplitudes ψ 1 and ψ 2 ; refer to  which are the known values of the macro-temperature ϑ on the top and the left edges of the composite (area 1 and area 2; refer to Figure 6) and the known values of the fluctuation amplitudes of the temperature ψ1 and ψ2 on all of the edges of the composite (areas 1 to 4; refer to Figure 7). Thus, a system of non-uniform equations was formulated. To solve this system, some variation of the Finite Difference Method was used-the so-called Crank-Nicolson method (the mixed method). The characteristic parameter for this method was assumed to be equal to a half. The Crank-Nicolson method guarantees the convergence of the solution regardless of the division of the structure.

Results
Using the created algorithm of the Finite Difference Method [23,25,26,28], calculations were conducted for the three structures described above. The two-dimensional, nonstationary heat-conduction problem of composites with base dimensions equal to L1 = L2 = 1 [m] was analyzed. During the analysis, the composite was assumed to consist of ten cells. In the case of the periodic structure, the dimension of the first sub-cell γ2 is taken initially as equal to 0.25; however, in the case of the structure with a functional gradation of properties, the function of gradation γ2(x2) equals x2/L2. On the other hand, in reference to the biperiodic structure, the volume share of the first sub-cell is originally equal to γ1 · γ2, where γ1 = γ2 = 0. 25    Then, the average temperature distribution was compared in a homogeneous structure with properties of the first sub-cell (blue line), then the second sub-cell (orange line); the periodic structure where γ 2 = 0.25 (grey line); the periodic structure where γ 2 = 0.5 (yellow line); the biperiodic structure where γ 1 = γ 2 = 0.25 (light-blue line); and the biperiodic structure where γ 1 = γ 2 = 0.5 (green line). This comparison of the selected x 2 -coordinate (x 2 = 0.5 [m]) is shown in Figure 14. For the readability of the chart, it is limited to nodes i from 5 to 20.           The values of the averaged temperature, shown in Figure 14, are summarized in Table  1. The lowest temperature values (ignoring the homogeneous structure) were obtained for the periodic structure where γ2 = 0.5, and the highest were obtained for the biperiodic   The values of the averaged temperature, shown in Figure 14, are summarized in Table  1. The lowest temperature values (ignoring the homogeneous structure) were obtained for the periodic structure where γ2 = 0.5, and the highest were obtained for the biperiodic The values of the averaged temperature, shown in Figure 14, are summarized in Table 1. The lowest temperature values (ignoring the homogeneous structure) were obtained for the periodic structure where γ 2 = 0.5, and the highest were obtained for the biperiodic structure where γ 1 = γ 2 = 0.25. The distinctions in values of the averaged temperature ϑ between these structures are in the order of 25.3-77.1% relative to a higher value.

Periodic and Biperiodic Structure
Then, the fluctuation amplitude ψ 2 distribution was compared (fluctuation amplitude ψ 1 does not occur in the case of the periodic structure) in the periodic structure where γ 2 = 0.25 (grey line); the periodic structure where γ 2 = 0.5 (yellow line); the biperiodic structure where γ 1 = γ 2 = 0.25 (light-blue line); the and biperiodic structure where γ 1 = γ 2 = 0.5 (green line). This comparison for the selected x 1 -coordinate (x 1 = 0.5 [m]) is shown in Figure 15. For the readability of the chart, it is limited to nodes j from 1 to 16.  Then, the fluctuation amplitude ψ2 distribution was compared (fluctuation amplitude ψ1 does not occur in the case of the periodic structure) in the periodic structure where γ2 = 0.25 (grey line); the periodic structure where γ2 = 0.5 (yellow line); the biperiodic structure where γ1 = γ2 = 0.25 (light-blue line); the and biperiodic structure where γ1 = γ2 = 0.5 (green line). This comparison for the selected x1-coordinate (x1 = 0.5 [m]) is shown in Figure  15. For the readability of the chart, it is limited to nodes j from 1 to 16.  The lowest fluctuation amplitude values (the absolute values) were obtained for the biperiodic structure where γ 1 = γ 2 = 0.25 and the highest were obtained for the periodic structure where γ 2 = 0.5. The distinctions in values of the fluctuation amplitude ψ 2 between analyzed structures are in the order of 13.65-72.21% relative to a higher value.

Periodic and Functionally Graded Structure
In this section, the approximated results of the Finite Difference Method are presented in the form of 3D-maps of the averaged temperature ϑ and fluctuation amplitudes ψ 1 and ψ 2 in reference to the functionally graded structure; refer to confer Figures 16-18.

Periodic and Functionally Graded Structure
In this section, the approximated results of the Finite Difference Method are presented in the form of 3D-maps of the averaged temperature ϑ and fluctuation amplitudes ψ1 and ψ2 in reference to the functionally graded structure; refer to confer Figures 16-18.

Periodic and Functionally Graded Structure
In this section, the approximated results of the Finite Difference Method are presented in the form of 3D-maps of the averaged temperature ϑ and fluctuation amplitudes ψ1 and ψ2 in reference to the functionally graded structure; refer to confer Figures 16-18.   Then, as in the previous section, the averaged temperature distribution was compared in a homogeneous structure with the properties of the first sub-cell (blue line) and of the second sub-cell (orange line); the periodic structure where γ 2 = 0.25 (grey line); the periodic structure where γ 2 = 0.5 (yellow line); and the structure with a functional gradation of properties where γ 2 = x 2 /L 2 (light-blue line). As before, the comparison for the selected x 2 -coordinate and for nodes i from 5 to 20 is shown in Figure 19. Then, as in the previous section, the averaged temperature distribution was compared in a homogeneous structure with the properties of the first sub-cell (blue line) and of the second sub-cell (orange line); the periodic structure where γ2 = 0.25 (grey line); the periodic structure where γ2 = 0.5 (yellow line); and the structure with a functional gradation of properties where γ2 = x2/L2 (light-blue line). As before, the comparison for the selected x2-coordinate and for nodes i from 5 to 20 is shown in Figure 19. The values of the averaged temperature, shown in Figure 19, are summarized in Table  2. For most of the nodes (ignoring the homogeneous structure), the highest values of the averaged temperature were observed for the functionally graded structure. The distinctions in values of the averaged temperature ϑ between the structure with a functional gradation of properties and the periodic structure where γ2 = 0.5 are in the order of 25.91-55.22% relative to a higher value.  Then, as in the previous section, the averaged temperature distribution was compared in a homogeneous structure with the properties of the first sub-cell (blue line) and of the second sub-cell (orange line); the periodic structure where γ2 = 0.25 (grey line); the periodic structure where γ2 = 0.5 (yellow line); and the structure with a functional gradation of properties where γ2 = x2/L2 (light-blue line). As before, the comparison for the selected x2-coordinate and for nodes i from 5 to 20 is shown in Figure 19. The values of the averaged temperature, shown in Figure 19, are summarized in Table  2. For most of the nodes (ignoring the homogeneous structure), the highest values of the averaged temperature were observed for the functionally graded structure. The distinctions in values of the averaged temperature ϑ between the structure with a functional gradation of properties and the periodic structure where γ2 = 0.5 are in the order of 25.91-55.22% relative to a higher value. The values of the averaged temperature, shown in Figure 19, are summarized in Table 2. For most of the nodes (ignoring the homogeneous structure), the highest values of the averaged temperature were observed for the functionally graded structure. The distinctions in values of the averaged temperature ϑ between the structure with a functional gradation of properties and the periodic structure where γ 2 = 0.5 are in the order of 25.91-55.22% relative to a higher value.
The fluctuation amplitude ψ 2 distribution was also compared in the periodic structure where γ 2 = 0.25 (grey line); the periodic structure where γ 2 = 0.5 (yellow line); and the functionally graded structure where γ 2 = x 2 /L 2 (light-blue line). This comparison for the selected x 1 -coordinate and nodes j from 1 to 16 is shown in Figure 20.  The fluctuation amplitude ψ2 distribution was also compared in the periodic structure where γ2 = 0.25 (grey line); the periodic structure where γ2 = 0.5 (yellow line); and the functionally graded structure where γ2 = x2/L2 (light-blue line). This comparison for the selected x1-coordinate and nodes j from 1 to 16 is shown in Figure 20.

Convergence Study
For each of the considered structures, a solution convergence analysis was carried out. The analysis consisted of examining the differences in the values of the unknowns. The solution was assumed to converge when, by increasing the number of intervals Δx1 The lowest fluctuation amplitude values (the absolute values) were obtained for the periodic structure where γ 2 = 0.25 and the highest were obtained for the periodic structure where γ 2 = 0.5. The distinctions in values of the fluctuation amplitude ψ 2 between analyzed structures are in the order of 34.88-69.41% relative to a higher value.

Convergence Study
For each of the considered structures, a solution convergence analysis was carried out. The analysis consisted of examining the differences in the values of the unknowns. The solution was assumed to converge when, by increasing the number of intervals ∆x 1 and ∆x 2 , the differences in the values of the unknowns decreased and were not significant. The results shown in Sections 5.1 and 5.2, obtained for the number of intervals ∆x 1 and ∆x 2 , were equal to 40.

Biperiodic Structure
At this point, the results of a solution convergence test for the biperiodic structure are shown. In Table 3, the values of the averaged temperature ϑ(x 1 ,x 2 ) and the fluctuation amplitudes ψ 1 (x 1 ,x 2 ) and ψ 2 (x 1 ,x 2 ) for a selected point (x 1 = L 1 /2, x 2 = L 2 /2) and changing number of intervals ∆x 1 and ∆x 2 are summarized. In addition, the differences in values of the averaged temperature with the smallest increase in the number of intervals were calculated. These differences were denoted as ∆ϑ(x 1 ,x 2 ). Analogous differences were calculated for fluctuation amplitudes and denoted as ∆ψ 1 (x 1 ,x 2 ) and ∆ψ 2 (x 1 ,x 2 ). As can be observed in Table 3, the differences in the values of the unknowns decrease as the number of intervals increases. The differences with the smallest increase in the number of intervals are in the order of 0.35-1.02% relative to a higher value.

Periodic Structure
At this point, the results of a solution convergence analysis for the periodic structure are shown. In Table 4, the values of the averaged temperature ϑ(x 1 ,x 2 ) and the fluctuation amplitude ψ 2 (x 1 ,x 2 ) for a selected point (x 1 = L 1 /2, x 2 = L 2 /2) and changing number of intervals ∆x 1 and ∆x 2 are summarized. As in the case of the biperiodic structure, the differences in values of the averaged temperature ∆ϑ(x 1 ,x 2 ) and the fluctuation amplitude ∆ψ 2 (x 1 ,x 2 ) were calculated and analyzed.
The differences in the values of the unknowns decrease as the number of intervals increases, and the smallest increase in the number of intervals are in the order of 0.45-1.06% relative to a higher value.

Functionally Graded Structure
For the functionally graded structure, the test of a solution convergence was carried out in the same way as in relation to the previous structures. In Table 5, the values of the averaged temperature ϑ(x 1 ,x 2 ) and the fluctuation amplitude ψ 2 (x 1 ,x 2 ) for a selected point (x 1 = L 1 /2, x 2 = L 2 /2) and changing number of intervals are summarized. The differences in values of the averaged temperature ∆ϑ(x 1 ,x 2 ) and the fluctuation amplitude ∆ψ 2 (x 1 ,x 2 ) were calculated and analyzed. The differences with a smallest increase in the number of intervals are in the order of 0.41-1.07% relative to a higher value.

Conclusions
The main goal of this work was to analyze three types of structures-periodic, biperiodic and those with a functional gradation of properties-under specific thermal boundary conditions. All these studies were possible thanks to a self-created algorithm, the Finite Difference Method, which was implemented into MAPLE software and applied to derived averaged model equations. This work shows the influence of the structure of the composite on the temperature distribution. The dimension of the sub-cell of the composite's cell is as important as the way and direction of placement of the components; therefore the volume shares of the materials with different properties are also important. As shown in Figures 14 and 19, the values of the averaged temperature ϑ vary significantly, with relative differences in the order of up to 77%. Similarly, as shown in Figures 15 and 20, the values of the fluctuation amplitude ψ 2 , which impacts in a major way on values of the total temperature field, vary substantially, with relative differences in the order of up to 72%. This leads to the conclusion that by designing the structure of the composite, the expected temperature distribution can be obtained, which in turn is important for the application of this type of construction. The analyzed structures also differ, importantly, in terms of self-weight or stiffness, which can be decisive in choosing the way and direction of the placement of the components and their volume shares.